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Abstract 

The last decade has seen an explosion in models that describe phenomena in systems medicine. Such 
models are especially useful for studying signaling pathways, such as the Wnt pathway. In this chapter 
we use the Wnt pathway to showcase current mathematical and statistical techniques that enable mod¬ 
elers to gain insight into (models of) gene regulation, and generate testable predictions. We introduce 
a range of modeling frameworks, but focus on ordinary differential equation (ODE) models since they 
remain the most widely used approach in systems biology and medicine and continue to offer great 
potential. We present methods for the analysis of a single model, comprising applications of standard 
dynamical systems approaches such as nondimensionalization, steady state, asymptotic and sensitivity 
analysis, and more recent statistical and algebraic approaches to compare models with data. We present 
parameter estimation and model comparison techniques, focusing on Bayesian analysis and coplanarity 
via algebraic geometry. Our intention is that this (non exhaustive) review may serve as a useful starting 
point for the analysis of models in systems medicine. 

Keywords: Wnt signaling, model development, nondimensionalization, asymptotic analysis, parameter 
inference, algebraic methods, model selection. 
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1 Introduction 

Despite the growing number of therapeutic options available to clinicians, gaps remain in our fundamen¬ 
tal understanding of many biological processes. Acquiring this additional knowledge requires that we 
focus on the molecular players that operate in intercellular and intracellular environments. Revealing the 
complex networks and dynamics that control cellular, tissue- and host-level behavior, may enable us to 
improve existing treatments and design new drug targets. 

Many intercellular signals are initiated by signaling proteins such as cytokines and hormones. When 
cytokines bind to receptors of a target cell, they trigger a cellular response by signal transduction path¬ 
ways: multistep sequences of intracellular signaling events and communication between molecules. 
Most of these molecules are proteins. Enzymes such as kinases and phosphatases, for example, cat¬ 
alyze (respectively) the addition/removal of a phosphate group to/from a substrate, and thus perform a 
crucial role in relaying information [1]. Phosphorylation (the addition of a phosphate group) can be as¬ 
sociated with protein activation, and information can be communicated downstream, engaging multiple 
signaling cascades by successive chemical reactions. While some reactions are linear, with the output 
proportional to the input [2], many are complex, involving feedback loops or pathway redundancies. 
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Often the output of these pathways is aetivation or inhibition of regulatory proteins ealled transeription 
faetors, whieh modify gene transeription and the eellular state. 

To turn a gene on, an aetivated transeription faetor transloeates from the eytoplasm into the nueleus, 
binds to the enhaneer or promoter region of DNA, and RNA polymerase transeribes the DNA template 
to synthesize RNA. Then messenger RNA (mRNA) leaves the nueleus and enters the eytoplasm where 
ribosomes translate mRNA into protein [1]. Conversely, transeription faetors may turn a gene off by 
repressing the reeruitment of RNA polymerase. These possible responses thus regulate protein synthesis. 
In addition to subeellular proeesses that ehanges in protein synthesis stimulate, proteins may be released 
by the eell and aet as signaling moleeules in other pathways. 

Gene regulatory pathways are erueial to the normal funetioning of eells, with many diseases eaused 
by dysfunetion of one or more pathways. For example, signaling pathways sueh as NF-fcB, MAP Kinase, 
and Wnt/j8-eatenin are involved in a host of eellular proeesses and funetions, ineluding eaneer. Due to 
their eomplexity, a systems approaeh is needed to understand normal and aberrant pathway funetion. 
Only by building theoretieal models that deseribe how eells signal and validating/updating them using 
experimental data ean we develop new drug therapies that target speeifie diseases. 

The remainder of the ehapter is organized as follows. In Seetion 2, we review methods used to model 
signal transduetion pathways, and introduee an exemplary enzyme kineties model. We then deseribe 
the biology of Wnt signaling, with referenee to relevant models, and introduee two models of the Wnt 
signaling pathway that we foeus on throughout the ehapter to demonstrate various teehniques. In Seetion 
3, we detail methods that ean be used to analyze a partieular model and diseuss the insight that eaeh 
approaeh ean generate. In Seetion 4, we introduee teehniques that ean be used to eompare models, 
ineluding some new methods for systems medieine. We eonelude in Seetion 5 with a diseussion of the 
different teehniques, and ideas for their further applieation in systems medieine. 

2 Mathematical Modeling 

Signaling pathways are eomplex and may be diffieult to understand by linear logie alone. Theoretieal 
models ean be used to gain insight into the dynamies of multiple bioehemieal interaetions. Construeting a 
mathematieal model is a nontrivial task, that requires suffieient understanding of the system to determine 
not only the type of model that should be used to address a partieular question but also the limitations of 
the model. After reviewing some of the modeling approaehes that are used to study signaling pathways, 
we foeus on ordinary differential equation (ODE) models. We introduee basie prineiples that ean be used 
to eonstruet ODE models and illustrate them by referenee to enzyme kineties and two models of the Wnt 
pathway. 
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2.1 Modeling Approaches for Systems Medicine 

Many processes associated with systems medicine in general, and signaling pathways in particular, can 
be modeled. These include: gene/protein abundances; gene/protein interactions; abundances of cellular 
species; the effects of cytokines, chemicals, drugs or other interventions on system or tissue-level phe¬ 
nomena. Modeling strategies for systems medicine can be classified as either deterministic or stochastic; 
we describe stochastic approaches briefly here, since the methods introduced in later sections are gener¬ 
ally only applicable to deterministic systems. 

Deterministic approaches describe systems for which, given full details of the model (parameter 
values and initial conditions), its time evolution can be determined exactly. This means that if a system is 
restarted multiple times from the same initial state it will always return to the same future states. Ordinary 
and partial differential equations (PDEs) are two examples [3]. PDEs with two or more independent 
variables (e.g., space and time) are more flexible than ODEs, but their simulation and analysis can be 
computationally expensive. Deterministic methods provide accurate descriptions of population-level 
behavior if the population sizes are large enough that the effects of random fluctuations can be neglected. 

Stochastic approaches describe systems whose temporal evolution has unpredictable elements due 
to randomness somewhere in the system. They are popular for modeling biological systems where ran¬ 
domness and heterogeneity abound, and should be used when population sizes are small enough that 
fluctuations cannot be ignored. In most cases, population averages will be recovered from a stochastic 
model when the abundances become large enough. One can construct stochastic models of protein dy¬ 
namics with stochastic differential equations [4] (i.e., ODEs with noise terms - often Gaussian - added). 
Such models can be used to study the dynamics of species that fluctuate about a well-defined mean value. 

Stochastic modeling can also be developed via agent-based approaches [5,6]. Here, individual agents 
act according to a set of rules. Eor example, within a given pathway, a protein could be phosphorylated 
or dephosphorylated with probabilities that depend on its environment. Such a framework treats protein 
species very differently to differential equation methods: each protein is viewed as an autonomous agent 
and population dynamics emerge in a “bottom up” manner. Whilst such methods may appeal to our 
intuition about protein heterogeneity, the approach is limited since analyses are often computationally 
expensive. As such, agent-based models should be used when population-averaged models fail to capture 
the behavior that the modeler seeks to describe. 

Cellular automata are a subset of agent-based models that impose spatial structure on the system by 
constraining the agents to lie on a grid, in two or three dimensions [7, 8]. The agents are updated via 
rules which may be deterministic or stochastic. Each grid point may be occupied by a finite number 
of cells (typically only one) and the model can accommodate multiple cell types. Cellular automata 
can account for spatial relationships between different cell types and have the advantage of being easy 
to interpret biologically. A challenge associated with these models is that the update rules may not 
translate clearly into biological hypotheses. Additionally, as for other agent-based models, simulation of 
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cellular automata can be computationally expensive. Fitting such models to data is at the limits of what is 
currently feasible since, despite significant advances in cellular imaging technology, obtaining cell data 
of sufficient resolution and quality to fit to a model is rare. 

The above overview of modeling approaches is not exhaustive: in limited space, we make no mention 
of Boolean, semi-quantitative, hybrid, or branching processes. Instead, we continue by explaining how 
to develop ODE models for signaling pathways. 


2.2 Formulating Mathematical Models of Signaling Pathways 

In this section our focus is on using ODEs to develop dynamic models of signaling pathways. Two basic 
principles are integral to the development of such models: 

• The Principle of Mass Balance states that the rate of change of a species is equal to the difference 
between the rate at which the species is added to the system and the rate at which it is removed; 

• The Law of Mass Action states that a reaction proceeds at a rate proportional to the product of its 
reactants. 

If, for example, substrate A is irreversibly phosphorylated by enzyme B, to produce C then we write 

A + B^C + B, (2.1) 


where ri is the rate at which phosphorylation occurs. We construct ODEs that describe the dynamics of 
A, B and C by appealing to the Principle of Mass Balance and the Eaw of Mass Action: 


dA 

dt 


-r\AB, 


dB 

dt 


—riAB + r^AB = 0 , 


dC 

dt 


riAB. 


( 2 . 2 ) 


By inspecting the above ODEs, it is straightforward to deduce that the following quantities are preserved: 


A C — Aq A Cq , and B — Bq , 

where A{t = 0) = Aq, B{t = 0) = Bq and C{t = 0) = Co are prescribed as initial conditions. We can 
exploit these Conservation Laws to simplify the governing equations: in this case, we can eliminate both 
B and C and our model reduces to give 

dA 

— =—riBoA, with A(t = 0) = Aq A(t) = 


Thus, substrate levels decay exponentially, at rate riB^. 



6 


Case Study I: The Enzyme Kinetics Model. 

We now consider a biochemical reaction that is catalyzed by an enzyme. In more detail, the enzyme 
E binds reversibly with the substrate S to form a complex C. While complexed with the substrate, the 
enzyme converts it into a product P and the enzyme is recovered. We represent these reactions as follows: 

E + S^c'^E + P. 

k-i 


By applying the Law of Mass Action to this reaction scheme and appealing to the Principle of Mass 
Balance, we deduce that the following system of ODES describe the time-evolution of S, E, C and P: 


dS 

dt 

dE 

dt 

dC 

dt 

dP 

dt 


-kiES + k-iC, 


-kiES + {k-i+k2)C, 


kiES-{k-i+k2)C, 


k2C. 


(2.3) 

(2.4) 

(2.5) 

(2.6) 


If we assume further that S{t = 0) = Sq, E{t = 0) = Eq, C{t = 0) =0 and P{t = 0) = 0, and take suitable 
combinations of the governing ODEs then we deduce 


d 

dt 


{E + C)=0 and 


d 

dt 


{S + C + P)=0, ^E + C 


Eo and S + C + P = So, 


We can exploit these conservation laws to eliminate E and P and obtain the following reduced model: 


dS 

— = -ki5(Eo-C)+k_iC, 

(2.7) 

— =ki5(Eo-C)-(k_i+k2)C, 
dt 

(2.8) 

with S{t = 0)=So and C(t = 0)=0. 

(2.9) 


2.3 Modeling Wnt Signaling 

Wnt signaling is implicated in many biological processes. The pathway is activated when Wnt ligands 
bind to specific receptors on the cell surface, resulting in the stabilization and nuclear accumulation of the 
transcriptional co-activator j3 -catenin. Canonical Wnt signaling encompasses cellular responses to exter¬ 
nal Wnt stimuli mediated by j3-catenin. Non-canonical describes cellular signaling and responses to Wnt 
not mediated by j3-catenin. The canonical Wnt pathway plays a key role in essential cellular processes 
ranging from proliferation and cell specification during development to adult stem cell maintenance and 
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Membrane 


Da ^ Di H Yap/Taz 



Figure 1. Reaction scheme that incorporates many different Wnt signaling models and additional 
molecular players (e.g., Yap/Taz). Solid arrows denote direct reactions; long-dashed arrows denote 
species that act as catalysts in degradation reactions; and dotted arrows denote alternative paths for the 
direct activation of Y. Note that active/inactive forms of Y are equivalent to active/inactive forms of 
ANG. 

wound repair [9]. Dysfunction of Wnt signaling is implicated in many pathological conditions, including 
degenerative diseases and cancer [10-12]. Despite further molecular advances [13-15], certain details 
of the dynamics of the pathway are still not well understood. 

The basic steps that constitute canonical Wnt signaling are as follows (although these are not undis¬ 
puted; discussed below): Wnt binds to cell-surface receptors Frizzled and LRP5/6 [11] that transduce a 
signal via a multi-step process involving Dishevelled (Dsh) to the so-called destruction complex (DC). 
The DC contains forms of Axin, adenomatous polyposis coli (APC), glycogen synthase kinase 3 (GSK- 
3), and casein kinase la (CKla). In the absence of a Wnt signal, the DC actively degrades j3-catenin 
- which is being continually synthesized in the cell - by binding and phosphorylating the protein and 
thus marking it for proteasomal degradation. Following Wnt stimulation, degradation of j3-catenin is 
inhibited through phosphorylation of DC member proteins. This leads to accumulation in the cytoplasm 


































of free j8-eatenin, which is able to translocate to the nucleus where it can form a complex with T-cell 
factor (TCF) and lymphoid-enhancing factor (LEF) proteins and, thereby, influence the transcription of 
target genes associated with processes such as self-renewal and proliferation [16,17]. 

In addition to these core mechanisms, evidence for other important processes has been found, some 
of which may challenge the Wnt signaling paradigm. Spatial localization within the cell has been found 
to be important not only for j3-catenin but also for Dsh and DC member proteins including Axin, APC, 
and GSK-3 [18-23]. There is also evidence of competitive binding of j3-catenin to cell membrane pro¬ 
teins such as E-cadherin [24] and intricate cross-talk with the Hippo pathway, this being mediated by 
Yap and Taz which promote translocation of cytoplasmic j3-catenin to the nucleus via phosphorylation 
and then compete with TCE for j3-catenin in the nucleus [25] . This spatial organization of Wnt pathway 
members may be key to understanding the pathway, as some modeling suggests [26,27]. Equally, an 
alternative description for the degradation of j3-catenin exists: in this picture, j8-catenin can be actively 
degraded while still bound to the DC, rather than being released marked for degradation [28]. Discrim¬ 
inating between competing hypotheses is needed in order to fully elucidate canonical Wnt signaling: 
mathematical modeling is a natural framework within which to achieve this. 

The first quantitative model of Wnt/j3-catenin signaling was developed in 2003 [29], based on data 
from Xenopus extracts. Eormulated as a system of DDEs, the model describes known interactions be¬ 
tween core components of the canonical pathway, these being Wnt, Dishevelled, GSK3j8, APC, Axin, 
j8-catenin and TCP. The DC is assumed to act only in the well-mixed cytoplasm and, hence, only cy¬ 
toplasmic levels of pathway components are considered. Since its publication, the Eee model has been 
extended in many ways (for recent reviews of mathematical models of Wnt signaling, see [16,30]). The 
effect of mutations in APC was investigated by Cho et al. [31], the action of Wnt inhibitors was studied 
by Kogan et al. [32], and the impact of Wnt-ERK cross-talk considered by Kim et al. [33]. The effect 
of competition for j3-catenin with adhesion proteins was investigated by van Eeeuwen et al. [34] while 
Schmitz showed how shuttling of core proteins between cytoplasm and nucleus could influence pathway 
dynamics [35,36]. More recently, a new shuttling model was constructed that accounts not only for ex¬ 
change of pathway proteins between the nucleus and cytoplasm, but also degradation of j3 -catenin while 
it is bound to active destruction complex (DC) and activation of the DC by dephosphorylation of its 
components [27]. Table 1 summarizes the key features of some of these models and Pigure 1 illustrates 
the localization and known interactions between key proteins involved in Wnt signaling. 

We now present the Eee model [29] and the Schmitz model [36], using the notation presented in Table 
2. These models, together with the enzyme kinetics model presented above, will be revisited throughout 
the chapter to illustrate how the techniques introduced in sections 3 and 4 are applied to specific models. 
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Biological Feature 

Lee 

van Leeuwen 

Schmitz 

Shuttle 

j3 -eatenin produetion 

/ 

/ 

/ 

/ 

j3 -eatenin degradation 
(independent of DC) 

/ 

/ 

X 

/ 

j3 -eatenin degradation (dependent 
on DC) 

/ 

/ 

/ 

/ 

j3 -eatenin sequestration by DC 

/ 

/ 

/ 

/ 

j3-eatenin sequestration by APC 

X 

X 

X 

X 

Shuttling of speeies between 
eytoplasmie and nuelear 
eompartments 

X 

X 

/ 

/ 

Aetivation/inaetivation of DC 

/ 

/ 

/ 

/ 

Interaetion with adhesion 

moleeules 

X 

/ 

X 

X 

Two j3-eatenin forms: transeription 
only and transeription or adhesion 

X 

/ 

X 

X 

DC is represented by its eonstituent 
parts 

/ 

X 

X 

X 

j8-eatenin binds individual parts 
APC and Axin as well as DC 

/ 

X 

X 

X 

j3 -eatenin binds to TCP to promote 
transeription of target genes 

/ 

/ 

/ 

/ 


Table 1. Comparison of features aeross different models of Wnt signaling. For further details 
see [27,29,34,36]. 


Case Study II: The Lee Model 

In its original form, the Lee model eomprises 15 time-dependent ODEs for protein speeies and eomplexes 
that partieipate in the Canonieal Wnt pathway, the reaetion rates being based on mass aetion kineties [29]. 
The model targets the assembly of the destruetion eomplex from the eonstituent parts of APC, Axin and 
GSK3j3. It does not distinguish between nuelear and eytoplasmie eompartments, instead assuming that 
all speeies are uniformly distributed throughout the eell. A sehematie diagram of the reaetions deseribed 
in the Lee model is given in Figure 2. Using the variable names defined in Table 2 and primes to denote 




10 


Symbol 

Species 

Forms 

X 

j3 -catenin 

Xp - marked for proteasomal degradation 

Y 

Destruction complex 

Ya - active 


(APC/Axin/GSK3j8) 

Yi - inactive 

D 

Dishevelled 

Da - active 

Di - inactive 

A 

APC 


N 

Axin 


G 

GSK3j8 


T 

TCP 


C 

Complex 

Cxy - complex of X and Y (etc.) 


Table 2. Definition of notation for the variables used by the Lee and Schmitz models. 


differentiation with respect to time, the DDEs that specify this model are: 

D/ = -ttiA + aiDa, (2.10) 

Dj = aiDi-a2Da, (2.11) 

yj = — OtaJa — CCloXYa + anCxY + OC^CxYp, (2.12) 

Yi = o^GCna - ocsDaYi - a^Yt + a4Ya - ajYt, (2.13) 

G' = a^DaYi — a^^GCxA + o;?!)-, (2.14) 

Cna = ocsDaYi — UfiGCxA + 0CiYi + a^NA — (X<)Cxa , (2.15) 

a' = — cc^NA + cxtjCxfA — 0t2i2fA + cc22CxA') (2.16) 

Cxy' = ocioXYa — auCxY — ciiiCxy, (2.17) 

CxYp = otnCxY — oti'iCxYp, (2.18) 

Xp = ccnCxYp — ciuXp, (2.19) 

X' = —axoXYa + otiiCxy + 0:15 — (Xi(,X — ai^XT + <X2oCxt — CI 21 XA + OC 22 CXA, (2.20) 

N' = -a^NA + ttgCwA + 0:17 - otisA^, (2.21) 

T' = —ttigXr + o:2oQvr, (2.22) 

CxT =ai9XT-a2oCxT, (2.23) 

CxA = a2xXA-a22CxA- (2.24) 


To facilitate comparison with the Schmitz model (see below), the non-negative rate constants a^., k G 
(1,2, ...,22) have been redefined from those used in [29]. Wnt dependence is incorporated via the pa¬ 
rameter a\ = «! (VT) that controls the activation of Dsh. 
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Inspection of Eqs. (2.10)-(2.24) reveals that there are four conservations laws: 


Dq — Di + Da , 

Go = G + Yi + Ya+CxY + CxYp, 

Ao = A + Yi + Ya + CxY + Cxyp + Cxa+Cna, 
To = T + CxT , 


the constants Do,Go,Ao and Tq denoting the (assumed constant) levels of Dishevelled, GSK3j3, APC 
and TCP initially present in the system. These conservation laws are consistent with experimental ob¬ 
servations which suggest that levels of these proteins do not fluctuate during Wnt signaling (i.e. they are 
produced and degraded at the same rates). They can be used to eliminate 4 variables and, in so doing, to 
reduce the model from 15 to 11 DDEs. Purther simplifications are achieved by assuming that all binding 
processes, except those for the binding of GSK3j3 to APC/Axin, reach equilibrium rapidly and that all 
species involving axin are present at low levels. Under these assumptions, and after some algebra, the 
following expressions for Dq,G,A, T,Xp,CxT,CxYp and Cna are obtained: 


Di = Do—Da, G = Go, A = 


Ao 


1 + 


T = 


To 


1 + 


ttl 9 V ' 

aio^ 


Xp = —CxY, 
ai4 


CxT = 


X.To 

l + f^X’ 
«20 


CxA = 


AqX 

1 + f^X’ 
«22 


CxYp = —CxY 

ai3 


^ _ ag AoN 


and a reduced system of 7 ODEs for the remaining species is eventually recovered (equations not pre¬ 
sented since they are rather involved and less instructive than Eqs. (2.10)-(2.24)). In [29], and [37] this 
model reduction is performed in an ad-hoc manner; it would be instructive to repeat it by first nondimen- 
sionalizing the governing equations (see section 3.2) and using asymptotic analysis to perform the model 
reduction (see section 3.3). 


Case Study III: The Schmitz Model 

Pike the Pee model, the Schmitz model [36] focuses on the canonical Wnt pathway. Key differences 
between the Pee and Schmitz models are that the latter distinguishes between the cytoplasm and nucleus 
and accounts for exchange of j3 -catenin and DC between these compartments (see Table 2 and Pigure 3 
for further description). In each compartment, DC binding to j3-catenin leads to its phosphorylation, and 
phosphorylated j3 -catenin is degraded. We use subscript n to denote species residing in the nucleus with 
the exception of TCP (T) and the jS-catenin-TCP complex (Cxt)', since these species are localized in the 
nucleus and to facilitate comparison with the Pee model, the subscript is omitted. Using notation that is 
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Figure 2. Schematic of the Lee model [29], which describes the activation of the destruction complex 
and its effect on j8-catenin in a single cellular compartment (cytoplasm and nucleus combined). 
Notation of the model species is given in Table 2. Solid arrows represent reactions and dashed arrows 
represent catalytic processes. 

modified from that used in [36], the ODEs that define fhe Schmifz model are: 


= ^ + (SzX^ - SiX) + (SeCxr - S 5 XY,), (2.25) 

X'„ = { 8 iX - & 2 Xn) + {d^CxYn - d^X„Yan) + {5nCxT - ^iiX^T), (2.26) 

X'p = &iCxY-5nXp, (2.27) 

X'p„ = 8 wCxYn-duXpn, (2.28) 

Y'a = {^AYan - W + (SeCxY - dsXYa) + &jCxY + - 5i5T,), (2.29) 

y; = SisYa - SieYi, (2.30) 

C = (53F« - 84Ya„) + (S^CxYn - 8 ^X„Yan) + SyoCxYn. (2.31) 

C'xY = {SsXYa - SeCxY) + SjCxy, (2.32) 

= (SsX„Ya„ - SgCxYn) - SioCxYn, (2.33) 

T'= SuCxr - SnX„T, (2.34) 

c;,r = duXnT-SnCxr, (2.35) 


where Sn (k = 1,2,... ,17) are non-negafive rafe consfanfs and 5i5 = 5i 5 (IT) so fhaf Wnf acfs fo inactivate 
fhe desfrucfion complex in fhe cyfoplasm. 
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Cytoplasm 

Xp^{D 



Figure 3. Schematic of the Schmitz model [36], which describes the interaction between j8-catenin and 
the destruction complex in two cellular compartments: cytoplasm and nucleus. Notation of the model 
species is given in Table 2. 


By taking appropriate combinations of Eqs. (2.25)-(2.35), it is straightforward to show that there are 
two conservation laws: 


F, + Ya + Yan + CxY + CxYn = YjoT and T + CxT = 7b, (2.36) 

the constants To and Tq denoting, respectively, the total number of molecules of DC and TCF in the 
system, as determined from the initial conditions. These identities may be used to reduce the order of 
the Schmitz model from 11 to 9. As explained below, further systematic simplifications may be possible 
following model nondimensionalization and parameter estimation. 

3 Techniques for the Analysis of a Specific Model 

Once model construction is complete, the modeler aims to extract from it new insight. This can be done 
in a number of ways: if no data are available, standard mathematical techniques can be used to increase 
understanding of the behavior of the model; however if data are available, then it may be possible to 
estimate model parameters. In this section we describe a number of techniques, some standard and 
others less so, that can be used to analyze models. We demonstrate these methods by reference to the 
models of enzyme kinetics and Wnt signaling introduced in Section 2. 
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3.1 Steady State Analysis 

Broadly speaking, the behavior of an ODE model ean be eategorized as either transient or steady state. 
The latter deseribes the behavior at large timeseales {t —)■ oo). For systems that reaeh single valued (i.e. 
not oseillating) steady states, we refer to the long time values that system variables take as the fixed 
points. Mueh theory exists for the analysis of fixed poinfs, whieh ean be helpful in eharaeferizing model 
behavior and prediefing fhe effeels of perfurbafions [38]. We eonfinue by ealeulafing fhe steady slates for 
Ihe enzyme kinelies model and fhe Sehmilz model (similar analysis ean be performed for fhe Lee model 
bul fhe resulling expressions are ralher involved and Iherefore omitted). 

Case Study I: The Enzyme Kinetics Model (Steady State) 

Setting ^ = 0 in Eqs. (2.3)-(2.6), we deduee that our model for enzyme kineties evolves to the following 
unique, steady state solution: 

5 = 0, E = Eo, C = 0 andE = 5o. 

Thus, as expeeted, the reaetion proeeeds until all of the substrate S has been eonverted to produet P. 
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and X solves a quadratic of the form 

0 = ^X^ + ^X + ^ (3.2) 

where the constant coefficients and are functions of the model parameters. For physically 

realistic solutions, we require > 0. Therefore, we conclude that this model has at most two steady 
states and at most one of them may be stable. 

As models increase in complexity, the algebra usually prohibits the construction of analytical expres¬ 
sions for the steady-state solutions. In the following sections we present other methods that can be used 
to generate insight in such situations. 

3.2 Nondimensionalization 

When a mathematical model is first developed, the independent and dependent variables typically repre¬ 
sent physical quantities (e.g., protein levels) which are measured in dimensional units (e.g., protein levels 
may be measured as the number of molecules per unit volume or the number of molecules per cell). The 
model may also contain parameters which relate to physical processes (e.g., reaction rates, Michaelis- 
Menten constants) and are also dimensional (e.g., rates may be measured per second, per hour or per 
day). Nondimensionalization involves recasting the model in terms of dimensionless (or unit-less) vari¬ 
ables. This process is instructive for several reasons. First, the number of model parameters is typically 
reduced. Second, the resulting dimensionless parameter groupings can provide useful information about 
the system’s behavior. Further, if estimates of these parameters can be obtained and then compared, it 
is possible to identify physical processes that dominate on a particular timescale and, thereby, rationale 
to simplify the governing equations. We illustrate these concepts by nondimensionalizing the enzyme 
kinetics and Schmitz models. 

Case Study I: The Enzyme Kinetics Model (Nondimensionalization) 

We introduce the dimensionless variables z,s,e,c and p where 

t = TT, S = Sos, E = EQe, C = Eqc, P = Sop. 

and the timescale T is specified below. It is natural to scale the complex C with Eq since the amount of 
complex that forms is limited by the amount of enzyme present. If the enzyme is working effectively (i.e., 
serving as an efficient catalyst), then the amount of product created will be comparable to the amount of 
substrate. Ttherefore, we scale P with So rather than Eq. 

There are several possible choices for the timescale T. Consider Eq. (2.3). Initially, when C = 0, the 
maximum rate of uptake of S is kiEo and similarly the initial rate of uptake of E is kiSo- The associated 
timescales are Ti = l/(ki£'o) and T 2 = l/(ki5'o). Since enzyme levels are typically much smaller than 
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substrate levels (i.e., Eq/Sq = £ ^ 1), it is elear that T2/T1 = Eq/Sq <C 1. We eonelude that Ti represents 
a long timeseale, assoeiated with substrate depletion, while T 2 represents a short timeseale, assoeiated 
with the initial rapid uptake of enzyme. 

Resealing on the longer timeseale, so that t = TiZ = xj{kiEo), Eqs. (2.7)-(2.9) transform to give 


ds 

dx 

dc 
e — 


— 5(1 — c) + KeC, 
5(1 -c) - K,„C, 


5(t = 0) = 1, c(t = 0) = 0, 


where £ = —, 


k-i 

kiSo 


and Km 


k-i +^2 
kiS'o 


(3.3) 

(3.4) 

(3.5) 

(3.6) 


Comparing Eqs. (2.7)-(2.9) and (3.3)-(3.6) we note that nondimensionalization has redueed the number 
of model parameters from five to three. We remark further that in Eq. (3.4), the initial eonditions supply 
dc{0)/ dx = I/e. Thus, if £ <C 1 then c will initially inerease very rapidly on the timeseale X. 


Case Study III: The Schmitz Model (Nondimensionalization) 

The proeedure for nondimensionalizing the Sehmitz model is identieal to that used for the enzyme ki- 
neties model. As the dimension of the system inereases, and more proeesses are ineluded, the number 
of ways to reseale the independent and dependent variables inereases rapidly. In sueh situations, it is 
important to eonsider whieh variables are expeeted to vary and over what timeseale: the answers to these 
questions should help to identify appropriate sealings. 

When studying Wnt signaling, inaetivation of the DC plays a key role in the system dynamies and 
therefore when we nondimensionalize the Sehmitz model time is resealed so that t = x/ 8 i$ is 
the timeseale for inaetivation of the DC). Variables relating to free j3-eatenin (i.e X,Xn,Xp,Xpn) are all 
resealed with B = 80 / 815 , the amount of jS-eatenin produeed during the typieal timeseale t. This sealing 
eliminates 80 from the dimensionless equations (see below). When ehoosing the sealings for variables 
involving DC and TCE, we aim to preserve eonservation laws. Aeeordingly, guided by Eqs. (2.36), we 
seale Ya,Yi,Yan,CxY and CxYn with Yjot, the total amount of DC in the system. Similarly, we seale T 
and CxT with Tjqt, the total amount of TCE in the system. Summarizing, we have 

{X,X„,Xp,Xpn) = B{x,Xn,Xp,Xpn), (!« ; ? ^an •) CxYiCxYn) — YTOT{ya,yi,yan ; ^xyt^xyn) 


(r,Cxr) = 7b(0,c,e), 

where x(t),v„(t), .. .,Cxe{x) are dimensionless variables. Under these sealings, the Sehmitz model gives 
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the following nondimensional system: 


X = 
/ 

Xn = 

x' = 


^pn 

1 / 
^yan 

1 / 

— = 


ft) 


xy 


1 / 


-e' 

V 

1 / 

-c,e 


1 + {52Xn - 5i;c) + {S^Cxy - SsXJa), 

{5ix - SlXn) + (S9Cxyn - SsXnyan) + (SnC^e “ SuXnd), 

^loCxyn S\4Xpfi, 

^{S4yan - ^ya) + (SeCxy - Ssxya) + ByCxy + “ Ja), 

— (ya-Siey,, 

CO 

^4yan) “h (^C-xyn S^Xfiyan) “h ^lOCxyti: 

{§5Xya SfyCxy^ ^Cxyi 

{5gXnyan ^C-xyn) ^lOC-xym 
{BuCxT - SllXnO), 

{BiiXnd - BnCxo), 


(3.7) 

(3.8) 

(3.9) 

(3.10) 

(3.11) 

(3.12) 

(3.13) 

(3.14) 

(3.15) 

(3.16) 

(3.17) 


where primes denote differentiation with respect to z and 5, (/ = 1,..., 16) are the following dimension¬ 
less parameters: 


. 5i 

Ol = 

Ol5 


02 = 

0l5 


5 S^Ytot 5 ^Yjot 

06 = --, 07 = 


Bn = 


5iiro 


.^15 


B\2 — 


(0 = 


5o ’ 
5i22b 

(^)/gl5) 

Ytot 


^ 

Ol5 


. 54 

04 = 

Ol5 


58 = 


S^Ytot 


, B9 = 


.>15 


s 5 i3 
5i3 = 

5i5 


5 5i4 

5i4 = 

5i5 


5,= 


S 9 YTOT 


5io = 


SsYtot 

SiqYtot 


80 

? 5i6 

5i6 — 

5i5 


and 


V = 


(5o/5i5) 

To 


(3.18) 

(3.19) 

(3.20) 

(3.21) 


3.3 Asymptotic Analysis 

In applied mathematics, if the (dimensionless) governing equations contain a small parameter, it is com¬ 
mon to assume that there is an asymptotic expansion for the solution, as a power series in the small 
parameter. As we demonstrate below, this technique can be used systematically to simplify a mathemat¬ 
ical model and, in so doing, provide useful information about the dynamics of its components. 
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Case Study I: The Enzyme Kinetics Model (Asymptotics) 


A key assumption of the enzyme kinetics model is that initial enzyme levels are much smaller than 
substrate levels. This assumption is represented in the dimensionless model equations via the small 
parameter e = Eq/Sq 1. We exploit this small parameter by seeking a solution to Eqs. (3.3)-(3.5) of 
the form 

s{z) so{r) + esi{T), c(t) ~ co(t)+ eci(T). (3.22) 

Substituting with Eq. (3.22) in the governing equations and equating to zero terms of 0{e"), we deduce 
that, at leading order, and cq satisfy 


dso , 

— = tq,co-5o(l - Co), 
at 

0 = 5o(l - Co) - ^mCO, 
5o(0) = l, co(0) = 0. 


(3.23) 

(3.24) 

(3.25) 


Thus the ODE for c reduces to an algebraic relation, giving co in terms of ^o, and an ODE for with 
the implicit solution 


K„,logsoiT:)+so{T) =A-kt, co = 


So 


(3.26) 


^mESo 

where A is a constant of integration. A problem arises when we attempt to impose the initial conditions: 
it is not possible simultaneously to satisfy both initial conditions. This is because the leading order 
problem is of lower order than the original one. 

In order to resolve this problem, we use matched asymptotic expansions. We recall that c varies 
rapidly near T = 0 and, hence, examine the system dynamics near T = 0 by switching to the short 
timescale T = xle. In terms of T, the model becomes 


ds ^ 

— =e{KeC-s{l-c), 

dc 

— =s{\-c)-KmC, 

5(0) = 1, c(0) = 0. 


(3.27) 

(3.28) 

(3.29) 


where 5(r) = s{t),c{T) =c(t). As before, we seek asymptotic expansions for 5 and c in terms of £ <C 1, 
of the form specified at Eq. (3.22). In this way, we obtain the following leading order solutions for sq{T) 
and co(r): 

s~o{T) = \, co{T) = -—-. (3.30) 

“T i^m 

The above approximate solution is accurate near r = 0 but not for T = 0(1), whereas Eq. (3.26) is 
accurate for T = 0(1) but not for T <C 1 . The method of matched asymptotics involves choosing the 



19 


constant of integration A to match Eqs. (3.26) and (3.30) [39]. By imposing the matching conditions 

lim(^o(T),co(T)) = lim {SoiT),CciT)), 

T^o r^-oo 

we deduce that A = 1. 

In practice, similar asymptotic analyses can be used to study ODE models of signaling pathways. 
As we have seen, such models may involve large numbers of variables and parameters, and estimates 
for many parameters may be lacking. In such cases, progress can be made by using order of magnitude 
estimates for certain processes. Eor example, in [29], the authors assume that all binding reactions are 
rapid, apart from the binding of GSK3j3 to APC/Axin. Under this/a^t kinetics assumption, the ODEs for 
the relevant species reduce to algebraic equations, in the same way that, for the enzyme kinetics model, 
on the longer timescale the ODE for the complex c reduces to an algebraic relation (see Eq. (3.24)). 

To the best of our knowledge, the Schmitz model has yet to be subject to such asymptotic analysis. 
Referring to Eqs. (3.7)-(3.17), and by analogy with the asymptotic analysis of the enzyme kinetics model 
presented above, we note that the dynamics of the system will be strongly influenced by the ratios ft) and 
V. Eor example, if typical levels of j3-catenin are much greater than levels of TCE and DC then we 
could construct approximate solutions to the Schmitz model in the limit for which v ^ 1 <C ft). Such an 
analysis of the Eee model was performed by [40]. Since the details are rather involved, we summarize 
the key points below, and refer the interested reader to [40] for further details. 

Case Study II: The Lee Model (Asymptotics) 

Numerical simulations of the Eee model generated using parameter estimates reported in [29] (see Eig- 
ure 5) suggest that the processes involved in the Wnt signaling pathway act over at least two different 
timescales. Eee et al.’s parameter estimates indicate that the basal rate at which j3-catenin is degraded is 
much smaller than the rate at which the DC becomes inactive. This discrepancy is exploited to define a 
small parameter, rj = aig/ais, which is the ratio of the rate at which jS-catenin undergoes natural decay 
to the rate at which the DC becomes inactive. The dimensionless parameters are then rescaled by multi¬ 
plying them by appropriate powers of Tj so that they are 0(1). By retaining terms of leading order, the 
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Figure 4. Series of schematics showing which components of the Lee model of Wnt signaling are 
active on the short (top), medium (middle) and long timescales. The active components on each 
timescale are highlighted with bold borders. Figure reproduced from [40], with permission. 
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log Time (non-dim) 



(APC'/Axin'/GSK3p) 



Figure 5. Series of figures showing how the Lee model responds to a Wnt stimulus (IF = 1) that is 
applied at 1 = 0 when the pathway is in equilibrium (IF = 0) at f = 0. Also shown is the asymptotic 
solution obtained by matching the short, medium and long time approximations to the Lee model. 
There is good agreement between the approximate and numerical solutions at all timescales. Key: 
numerical simulations of the (dimensionless) Lee model, Eqs. (2.10)-(2.24) (solid line); Short, medium 
and long time approximations are represented by dash-dotted, dotted and dashed lines respectively. 
Figure reproduced from [40], with permission. 
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following reduced model is obtained: 

dDa , , 

-^ = a,W{\-Da)-a2Da, 

dF, ^ Ok,N 

— - -{a,Da + a2 + ai)Yi + Ya + j^-^, 

^d^ = aioXYa - ttiiCxF, 
at 

dY, 

dt I+K3X 

IdX _ _ 

- ^ = ai5 - aioXYa - aieX. 

T] di 


(3.31) 

(3.32) 

(3.33) 

(3.34) 

(3.35) 

(3.36) 


We remark that Eq. (3.31) decouples and if a constant Wnt stimulus is applied (W{t) = W, constant) 
then 

«! + a 2 

We note further that the time derivatives in Eqs. (3.31)-(3.36) are premultiplied by three different powers 
of 77 . This suggests that model processes act on three distinct timescales, a prediction that is consistent 
with the rapid fluctuations and slow increases depicted in Eigure 5. 

As for the enzyme kinetics model (Eq. 3.22), it is possible to analyze the reduced Eee model on the 
short, medium and long timescales for which t = 0(7]),0(1) and 0(r]^^) respectively. In each case, 
asymptotic expansions in powers of the small parameter 7 ] are sought and used to simplify the governing 
equations. The results of this analysis can be summarized as follows (see [40] for details). 


1. Short timescale (7 = 0(7])): all model variables except Y and Cxr are constant, at leading order. 
The dominant reaction is phosphorylation of j3-catenin by active destruction complex. 

2. Intermediate timescale (t = 0(1)): the dominant reaction is found to involve inactivation of the 
destruction complex. 

3. Eong timescale (t = 0(7]^^)): the dynamics are dominated by degradation of free j3-catenin. 


Pathway components acting on the short, intermediate and long timescales are highlighted in Eigure 4 
while Eigure 5 shows good agreement between the approximate solutions and those of the full model. 


3.4 Parameter Analyses 

The selection of model parameters, their physical meaning and numerical value are especially impor¬ 
tant; therefore, parameter analysis examines the response of the system to changes in parameters. Many 
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methods for estimating parameters depend on time eourse data. These data generally give a quantitative 
measure of the variable level, sueh as mRNA or protein eoneentration level, at different time points. 
Testing a model against experimental data is a good way to validate or invalidate it; however, gathering 
experimental data is often too expensive to determine all parameter values and overfitting, i.e., deseribing 
noise instead of the relationship, is a risk, as demonstrated for Wnt signaling later in this seetion. Fol¬ 
lowing parameter estimation (using optimization) or parameter inferenee (using statisties), a good way 
to test a model is by performing parameter sensitivity analysis: this evaluates qualitative or quantitative 
relationships between parameters and their effeet on the system outeome [41]. 

Parameter Estimation and Wnt Data 

Ultimately, every model should be tested against data, a proeess that ean either invalidate the model or 
provide evidenee in its favor, if it provides a good fit under aeeeptable eonditions. The aim is to esti¬ 
mate parameters that drive the model elose to the data; this ean be done using minimization teehniques. 
Effeetively, one ealeulates an objeetive funetion whieh is defined as fhe differenee befween fhe model 
simulated for parfieular value of paramefers K and fhe observafions (dafa), and aims fo minimize fhe error 
of fhe objeefive funetion, oflen performed iferafively [42^14]. 

Sinee publieafion of fhe Lee model [29], where esfimafes of fhe parameters eonfrolling Wnf signaling 
were based on dafa from Xenopus exfraefs, few sfudies have quanfifafively sfudied fhe dynamies of fhe 
Wnf pafhway. This knowledge gap means fhaf eurrenfly if remains diffieull fo fesf fhe models fhaf have 
arisen in reeenf years. This problem is nol uncommon in sysfems medicine. We also remark fhaf fhe 
Xenopus dafa gafhered by Lee et al. may be markedly differenl from fhose for mammalian Wnf signaling. 
In [13], dynamic changes in j3-cafenin levels were invesfigafed in Xenopus exfraefs. They demonsfrafed 
fhaf absolute levels of j3-catenin did nol diclale fhe Wnf signaling outeome: rafher fhe j3-calenin fold- 
change was fhe crucial variable. They used fhe Lee model fo fesf fheir experimenfal resulfs and, via 
sensifivify analysis, idenlified fhaf fhe model confirmed fheir experimenfal findings. 

Quanlificalion of Wnf signaling in mammalian cell lines was undertaken by [14,15]. Discrepancies 
wifh dafa from Xenopus exfraefs (such as higher Axin levels and lower APC levels in mammalian cells) 
highlighf fhe need for caution in dafa gafhering and for furlher quantification of fhe pafhway. Since fhese 
measuremenls were made al steady slale, Ihey do nol yel permil elucidalion of Iransienl Wnf signaling. 
More reeenf measuremenls of cytoplasmic and nuclear jS-calenin in response to a Wnf stimulus provide 
a valuable firsl look af fhe dynamics of fhe pafhway [45]. 

The above sfudies provide preliminary insighl into fhe Wnf pafhway bul much remains fo be done. 
The dafa are nol yel of sufficienl qualify to discriminate befween mosf models (which typically conlain 
many molecular species). Caulion musl be laken when applying dafa. Lor example, where dafa generated 
from non-mammalian sysfems may be used in a model fhaf addresses clinical oulcomes. Lor sysfems 
medicine fo have fhe grealesl impacl, modeling (wifh prediclion) and experimenlalion (fo fesf predictions) 
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must proceed iteratively. 

Parameter Inference 

There are often cases where it is either infeasible or impossible experimentally to determine values 
for parameters that describe a given model. In such cases, we may be able to estimate (some of) the 
parameters using statistical inference. In general the aim is to identify the values of the parameters, 6 , 
(ideally including corresponding confidence regions) for which a model best explains the data. 

A reliable way of doing so is to focus on the likelihood L{d), which is defined as fhe probabilify of 
observing fhe dafa (x) given paramefers (0): 


L{d) :=P{x\d). 

Varying d fo idenfify fhe value for which fhis probabilify is maximized gives fhe maximum-likelihood 
esfimafe. There is a rich liferafure on fhis topic and how confidence of fhe estimates can be assessed [46]. 

Likelihood esfimafes cenfer around fhe available dafa. In many circumsfances we may have addi¬ 
tional information, for example based on biophysical argumenfs, abouf which paramefer values can be 
ruled ouf. Incorporafing such prior information is hard in a pure likelihood framework, buf lies af fhe 
hear! of Bayesian inference [47]. Here inferences are based on fhe posterior distribution over model 
paramefers. The posterior disfribufion can be described sfarfing from Bayes rule: 

P{d\x)ocP{x\d)n{d). (331) 

P{Q\x), fhe probabilify of 0 given v, is called fhe posterior probabilify, P{x\Q) is fhe likelihood function, 
and 7r(0) is fhe prior probabilify (knowledge abouf paramefers before we begin tiffing fo dafa) [48]. As 
well as fhe full (join!) posferior disfribufion, one may also analyze fhe marginal posferior disfribufions 
which are fhe individual disfribufions over each paramefer. 

In cerfain cases, such as for large, complex sysfems, compufing fhe likelihood is impracfical. In such 
cases approximafe Bayesian compufafion (ABC) should be considered [49]. Insfead of fhe likelihood, a 
disfance funclion is used fo compare fhe acfual dafa wifh dafa simulafed by a model, denofed Xm- If the 
underlying model is given by / = /(x^l 0), then we express the ABC posterior function by 

Pabc{Q\x) l{A{x,Xni) < e)f{xm\d)n{d) (3.38) 

where A{a,b) denotes a distance measure between a and b and e is the tolerance level that determines 
how well real and simulated data should agree. 

By evaluating the posterior function, ABC allows the modeler to identify parameter regions that are 
of interest, and ignore those that are not. Furthermore, the posterior distribution gives information about 
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Figure 6. Data published in [45] were used to fit the Lee and Schmitz models using approximate 
Bayesian computation for parameter inference. j3-catenin concentration units were normalized based on 
their initial values. From the inference, we can see that the Lee model provides a better fit to the data. 


joint distributions in parameter space and can reveal multivariate dependencies between parameters. 

ABC for parameter inference has been implemented in the software package ABC-SysBio with sup¬ 
port for parallelization [50]. For the examples given below, we used the CUDA implementation of 
ABC-SysBio with a Euclidean distance measure between model and data [51,52]. Proceeding to analyze 
the Lee and Schmitz models, we do not try to infer ah of the model parameters, since this is not possible 
with the data available, but instead study a 3D subset of parameter space. We choose free parameters that 
have direct (or strong) influence on fhe dynamics of j3-cafenin, since fhis is fhe species for which we have 
experimenfal measuremenfs. The dafa used for hffing are published in [45]: fhey describe how fhe level 
of j3-calenin changes over lime in fhe cyloplasm and nucleus, following applicafion of a Wnl sfimulus fo 
fhe system. These dafa, alongside fhe resulfs of fhe paramefer inference, are shown in Figure 6. 

For fhe Lee model, we sfudy fhe jS-cafenin-DC binding rate (otio), thaf has a prior of [0,100], fhe 
j8-cafenin degradafion rale lhaf is independenf of fhe DC (otig), and fhe binding rale of jS-cafenin lo 
TCF (ai 9 ). The lafler Iwo paramefers bolh have priors of [0, Ij. The marginal poslerior dislribulions for 
Ihese Ihree parameters (Figure 7) show lhaf fhe j3-calenin-DC binding paramefer lakes values over fhe 
lower half of ifs prior range, whereas fhe ofher Iwo paramefers can lake any values spanning fhe prior 
range. This suggesls lhaf for fhis model fhe paramefer lhaf has fhe grealesl impacf on oufcome is fhe 
j3-calenin-DC binding rate, however we nole fhe larger prior range over fhis paramefer. 

For fhe Schmilz model, we sfudy fhe j3-catenin producfion rate (5o), the j3-catenin shuttling rate (5i), 
and the binding rate of j8-catenin to TCF (5ii). The prior used for each parameter is [0,1] and we see 
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Figure 7. Posterior distributions and sensitivity analysis for the Lee and Schmitz models. Histograms of 
marginal posteriors for each free parameter in the two models are shown. The marginal posterior is the 
probability distribution for a single parameter, given data describing j3 -catenin dynamics in cytoplasmic 
and nuclear compartments [45]. Principal component (PC) analysis allows us to assess the sensitivity of 
the parameters to small perturbations: the last PC, PC3, contains the most sensitive parameters. We see 
that for each model, two parameters dominate PC3 and, thus, are most sensitive in this system. 
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from Figure 7 that the marginal posterior distributions are relatively stiff: eaeh parameter is eonstrained 
to lie within a narrow range relative to its prior. In order to fit the data, the rates of jS -eatenin shuttling 
and binding to TCF must be low, whilst the rate of -eatenin produetion must be high. 


Sensitivity Analysis 

Sensitivity analysis investigates how a model responds to perturbations around a set of parameter values 
and eharaeterizes its robustness: a robust system is one for whieh perturbations of the parameters or initial 
eonditions do not ehange the outeome. However, many trade-offs between sensitivity and robustness 
exist [53-55]. 

Local sensitivity analysis determines how parameter perturbations affect the output of a system. 
Estimated or inferred parameters can be used as a baseline for parameter sensitivity. If the output of 
dx/dt = f{x,K) is approximated by a first-order Taylor series in a neighborhood of reference input 
values, then the local sensitivity coefficient, Sij is the partial derivative of the state to the parameter: 




dxi{t) 


(3.39) 


The elements si^j define a sensifivify mafrix S = dx/dK. This local mefhod provides informalion abouf 
fhe sensifivify in a given paramefer region buf nof fhe global sensifivify landscape. Local sensifivify 
analysis can reveal paramefers fhaf are sensitive or robusf fo perfurbafions in fhe region of inferesf. 

Principal componenf analysis (PCA) offers anofher way fo invesfigafe sysfem sensifivify. This fech- 
nique can be readily applied fo fhe posterior disfribufion obfained following Bayesian inference. The 
principal componenfs are consfrucfed by evaluafing fhe eigenvalues and eigenvecfors of fhe covariance 
mafrix of fhe paramefers: fhe firsf principal componenf (given by fhe largesf eigenvalue) corresponds fo 
fhe direcfion in which fhe posferior is mosf wide; fhe lasf principal componenf (given by fhe smallesf 
eigenvalue) corresponds fo fhe direcfion in which fhe posferior is mosf narrow [49, 56]. The lasf few 
principal componenfs represenf fhe mosf sensifive (or “sfiff” paramefers) [57]. 

In Figure 7, sensifivify analysis via PCA for fhe Lee and Schmifz models is shown. The principal 
componenfs (PC) are ordered 1-3 fhus PC3 is fhe lasf componenf and confains fhe mosf sensifive 
paramefer combinafions. Lor bofh models, PC3 is dominafed by fwo paramefers: fhe rafes of jS-cafenin 
binding fo fhe desfrucfion complex (DC) or fo TCP for fhe Lee model (aiOjOtig); and fhe rafes of j3- 
cafenin producfion or binding fo TCP for fhe Schmifz model (5o, 5ii). These resulfs suggesf fhaf fhe Lee 
model is more robusf fo changes in fhe j3-eatenin degradation rate (otie), and fhaf fhe Schmifz model is 
more robusf fo changes in fhe j8-eatenin shuffling rate (5i). 
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4 Techniques for the Comparison and Discrimination of Models 

Given a set of models that deseribe similar biologieal phenomena, a ehallenge is to determine whieh 
model best deseribes the system, given the evidenee available. In this seetion we deseribe two methods 
that enable eomparison and diserimination between models. The first employs approximate Bayesian 
eomputation, introdueed above, and has already gained a strong foothold in systems medieine [50,58- 
60]. The seeond is model diserimination with the use of algebraie matroids; as far as we know this is a 
reeent addition to the modeler’s toolkit and holds great potential for advanees in systems medieine. 

4.1 Model Selection via Approximate Bayesian Computation 

Returning now to the Lee and Sehmitz models, we eonsider how to ehoose between models using ABC 
model seleetion. We have already demonstrated how methods for parameter inferenee, sueh as ABC, ean 
yield the posterior distributions over the parameters of a model (given data) and diseussed briefly how 
this ean be interpreted. For two or more models (M, , / = !,...,«) some measure of the evidenee for eaeh 
model is needed [61], 

P{Mi\x)^P{x\Mi)n{Mi), (4.1) 

where (as previously) v represents the data, and n the prior probability. 

The ABC approaeh may be extended to parameter inferenee and model seleetion simultaneously 
using a joint spaee approaeh [49]. This may be performed for M models where M = 
by assigning to eaeh model (and parameters therein) a prior distribution and perturbation kernel that 
designates weights for model transition. The algorithm aeeepts N partieles at the Ef toleranee, whieh 
forms the joint posterior distribution P{a,M\x) and upon marginalizing over parameters, the marginal 
posterior distribution P{M\x) is approximated, providing a measurement for model seleetion. Bayesian 
model seleetion, like other approaehes ineluding the likelihood ratio test or Akaike Information Criteria 
(AIC), also penalizes over-parameterization. 

The AIC for model M,-, with / G {1,..., n}, is defined as 

AIC,' = -2logL{e*;x,Mi) + 2ki, (4.2) 

where L is the likelihood, and Q* and k, are (respeetively) the maximum likelihood parameter and num¬ 
ber of parameters in model M,. This eriterion, probably the best known model seleetion tool, makes 
explieit the penalty for an inereased number of parameters. However, as the amount of data inereases, 
the AIC introduees bias and tends to favor models that are over-parameterized. Therefore the Bayesian 
information eriterion (BIC), 


BIC; = -21ogL(e*;x,M;) +ki\ogn, 


(4.3) 
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Figure 8. Model selection via ABC for the Lee and Schmitz models. The results show that, over 
successive populations, evidence in favor of the Lee model grows until there is a high probability that 
this model will be selected, given the data published in [45]. 


may be preferred, as it remains unbiased for large samples, n. The BIC is effectively an approximation 
to the model probability (4.1); the penalty term, explicit in the AIC and BIC definitions, is implicit in 
(4.1), where it enters via the parameter priors for each model. 

Model selection chooses, from among a set of candidate models, the model that best explains ob¬ 
served data. Two things need to be kept in mind: (i) one model will always be chosen as the best but this 
does not mean that the model is necessarily a good one; ideally model selection should go hand-in-hand 
with model checking (and topological sensitivity analysis [62]). (ii) model selection depends on the data 
available for testing the different models; since different data may favor different models, careful exper¬ 
imental design should precede model selection. With these issues in mind we have the pragmatic choice 
about which statistical model selection framework to employ. Fully Bayesian, even in an ABC context, 
is more expensive than identifying the maximum likelihood parameter set and applying AIC or BIC. 

Shown in Figure 8 are the results of ABC model selection for the Lee and Schmitz models, with the 
probability of the model given for successive iterations (populations). We see that initially both models 
are equally probable, but subsequently the probability of selecting the Schmitz model drops to close to 
zero and we conclude that the Lee model is favorable given these data and parameter combinations. 








30 


4.2 Model Discrimination using Parameter-Free (Algebraic) Approaches 

When parameter values are unknown or cannot be estimated from data, one may still be able to discrimi¬ 
nate between competing models. We present two approaches, one that requires no data (rather qualitative 
insight into whether the system can have multiple responses) and another method which requires either 
highly resolved single cell data or multiple replicates of steady-state measurements. 

Precluding/Asserting Behaviors via Chemical Reaction Network Theory 

Chemical reaction network theory (CRNT) studies the structure of a model (which can also be described 
as a network) constructed from chemical reactions without relying on specific parameter values. The 
aim here is to use such theory to preclude (and sometimes assert) possible qualitative behaviors in the 
positive orthant, i.e., M>o. Cases where multiple positive states are stable (i.e., biologically accessible) 
are of particular biological importance for cellular decision making, for example, differentiation into one 
of two or more specialized cell lineages. 

The field of CRNT inifially focused on a sfrucfural properly of a model called deficiency, which could 
preclude mulliple sleady-slales [63, 64]. Then Iheorems were proved for precluding/asserfing multiple 
equilibria by sludying fhe cycles in Ihe graph of a nelwork, or fhe sign of fhe delerminanl of Ihe Jacobian; 
some of fhese approaches can provide condilions on fhe paramefers for behaviors such as bislabilily and 
oscillafions [65-70]. An excellenl and comprehensive survey of fechniques for mullisfafionarily was 
wrilfen by Joshi and Shiu [71]. One main fool for precluding mullisfafionarily of a model is lesling 
whelher if is injective (a model, including conservation relations, is injective if F{x,k) = F{x,k) 

X = x). Here we demonslrale fhe applicalion of mullisfafionarily lesls (developed for chemical reaction 
nefworks) lo Wnl signaling models. 

We begin wilh fhe Lee model. Firsl we lesl injeclivify, noting lhal while injeclivify precludes mul- 
lislalionarily, failure of injeclivify does nol imply mullisfafionarily. We use fhe algorifhms in fhe CRNT 
Toolbox lo determine whelher fhe system can ever admil multiple posilive sleady slales-mulfisfalionarify 
[72]. The Lee model fails injeclivify, bul cannol admil mulliple positive steady slates for any values of Ihe 
system parameters and/or lolal concenlralion amounls (algorilhms wilhin [72]). Conversely, Ihe Schmilz 
model has Ihe capacity for multiple sleady-slales; however, as calculated earlier, only one can ever be 
slable. Therefore, in Ibis example, since bolh models only can have one slable sleady-slale, il is diflicull 
lo use only qualilalive dala lo discriminate belween Ihem. Clearly, if dala suggested Iwo slable slates 
could exisl (for example via flow cylomelry), and all of Ihe dala had Ihe same initial conditions, Ihen one 
could rule bolh models out 





31 


Model Discrimination using Coplanarity via Algebraic Geometry 

When data from a model elearly supports a speeifie behavior — whether monostable, bistable, or os- 
eillatory, qualitative approaehes sueh as those mentioned above may be a good first step for elassifying 
models, espeeially if the data are not suffieient to estimate parameters. However if steady state data are 
available, then determining steady-state invariants may be helpful for determining whether a model is 
eompatible with given data using a statistieal parameter-free model diserimination method. 

Sinee often data are not available for all model speeies, variables must be eliminated. A system- 
atie teehnique from algebraie geometry proeeeds by eomputing the Grobner Bases of the model variety 
(studying the model at steady-state) and eliminating unobservable variables. The resulting steady-state 
invariant enables us to foeus on part of the system and to test whether the data suggests that the rela¬ 
tionships between speeies still hold. Notions of dependenee and independenee between model variables 
ean also be studied using algebraie matroids and were reeently applied to steady-state model diserimina¬ 
tion [27]. 

For smaller models, the steady states ean be determined explieitly. For example, for the Sehmitz 
model, the steady state values ean be expressed in terms of X and all other variables ean be eliminated 
by exploiting eonservation laws and using variable substitution (see Eqs. (3.1)-(3.2)). Either by hand, 
by eomputing the matroid, or by using Grobner bases, the polynomial relationship/algebraie dependenee 
between X and Xn in the Sehmitz model gives the following invariant: 

^ = ^>535456(58 -b 5i))X^ + (505257^9(55 -1-56) — 5i 535456(58 -I-59))AA„ — 5i 525/59(55 + 5^)X^, 


whieh vanishes at steady-state (i.e., = 0). Effeetively, we aim to test whether the data are eoplanar 

with our model, via the steady-state invariant transformation. Model eompatibility is determined by 
eomputing the eoplanarity error (A) via the singular value deeomposition of the matrix 


( 

^2 A2 XX„ 

V 


(h:\ 

hi = 0, 

\h3 ) 


where X denotes the observed value of speeies X. The null hypothesis (that the model is eompatible with 
the data) ean be rejeeted when the eoplanarity error (normalized smallest singular value) is less than a 
statistieal bound, whieh is determined by the Gaussian measurement noise in the data and the invariant 
strueture [73]. This method was reeently applied to j3-eatenin loealization data (eytoplasmic, A; and 
nuelear, X^) published in [27,45]. The Sehmitz model eould be ruled out if data were perturbed less than 
10^^ by measurement error/noise; for higher levels of noise, the model is eompatible. 
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5 Discussion 

Paradoxically, technological advances can sometimes create new challenges for clinicians. For example, 
as the number and variety of treatments for cancer increase, it can be difficult to identify the combination 
of treatments that will most benefit a given patient (if a unique, optimal treatment even exists). The 
situation is further complicated when we consider the different types of data that can be used as a basis 
for diagnosis and treatment planning; it is often impossible to integrate the available data by linear 
thinking alone. Systems medicine aims to address these challenges by developing mathematical and 
computational tools that integrate different types of information in order to generate objective decisions 
for patient treatment. In this chapter we have focused on ODE models, a class of models widely used in 
systems medicine, particularly to study signaling pathways. We have reviewed a variety of techniques 
that can be used to develop and analyze ODE models, using models of enzyme kinetics and the Wnt 
signaling pathway as test cases. 

Many of the techniques that we have presented are already well-established (such as model develop¬ 
ment, nondimensionalization, identification of steady state solutions, asymptotic analysis, and parameter 
sensitivity analysis); however others are less well-known (such as approximate Bayesian computation, 
chemical reaction network theory, and matroid-informed coplanarity). In addition to the benefit that these 
methods bring to the field, model developmenf for sysfems medicine - in ifs increasing sophisficafion - 
is helping fo sfimulafe furlher developmenf and applicafion of mafhemafical and sfafisfical fechniques. 

Many of fhe challenges in sysfems medicine arise because mosf biological processes, including pafh- 
ways, do nol acf in isolafion. Eor example, af fhe subcellular level, pafhway cross-falk can have a sig- 
nificanl effecf on cell funcfion. In particular, fhere is growing evidence of cross-falk befween Wnf and 
E-cadherin [74], Wnf and Erk [33]) and Wnf and fhe Hippo pafhway [75]. Even simplisfic models of 
such pafhway cross-falk quickly become large and demand sophisficafed techniques for fheir analysis. 
The sifuafion becomes more complex when we consider fhe impacf of signaling pafhways af fhe mul¬ 
ticellular and tissue scales. The impacf of Wnf signaling af fhe mulficellular and fissue levels has been 
sfudied Iheorefically, mosf prominenfly in models of infesfinal crypfs [76-79]. These models (for ex¬ 
ample) infroduce spatial dependence by imposing a graded Wnf disfribufion along fhe crypf axis [78] 
or provide comparison of a confinuum model wifh a cell-based model fhaf incorporafes heferogeneify 
and noise [79]. In [74], a multiscale model of inferacfions befween fhe pafhways affecfing j3-cafenin 
and E-cadherin is developed and used fo sfudy fhe role of epifhelial-mesenchymal fransifions in can¬ 
cer growfh and mefasfasis, whereas in [80] a simple rule-based model for cross-falk befween fhe Wnf 
and delfa-nofch pafhways is embedded wifhin discrefe epifhelial cell agenfs and used fo sfudy cell fate 
specification wifhin fhe infesfinal crypf. In addition fo fhese fheorefical sfudies (ever growing in com¬ 
plexify), more sophisficafed dafa collection is urgenfly needed as a basis for hypofhesis fesfing and model 
(in)validalion. 
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We end by proposing two grand challenges, whose solutions will bear much fruit in systems medicine. 
The first is to incorporate multiple levels of information — from biochemical reactions within a single 
cell, to tissue-level processes — into cohesive models. The second is to incorporate data which is re¬ 
solved in space and time into a theoretical framework. There are, of course, many other examples, and 
work in these areas should provide many exciting opportunities for theoreticians in systems medicine for 
years to come. 
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